This notebook shows a sample workflow for running hydrology simulations using the GSSHA rectangularly gridded simulator, supported by a suite of primarily open-source Python libraries. The workflow consists of:

Selecting parameters using widgets in a Jupyter notebook to control the model to simulate, including a watershed shape file.
Visualizing the watershed shape in a geographic context (projected into a suitable coordinate system and overlaid on map tiles from a web tile server).
If necessary, editing that watershed shape by hand and creating a new shape file with the edited result.
Selecting parameters to control the simulation, potentially overriding some selected earlier for the model creation (e.g. if running numerous conditions as a parameter sweep).
Visualizing and reviewing the inputs to the simulation.
Running the underlying simulation, collecting data on flood depth at each time point as well as the overall maximum flood depth per grid cell.
Visualizing the flood depth over time and the maximum flood depth.
Analyzing the simulation speed to help shape expectations for computational requirements for future runs.

Each of these steps is configured directly in this notebook, and can thus easily be scripted or iterated as needed. The set of parameters and precisely how they are configured is still being improved, and it can likely be made into a better match to users' needs in this domain. This workflow relies on fast raster regridding added to Datashader and exposed via HoloViews as part of the EarthSim project.

The underlying environment needed to run this workflow is set up as described in the README , and though already functional will need to be greatly simplified to be more usable and maintainable in practice. The workflow currently relies on downloading data from external servers that can be slow to access from some parts of the internet, so you may see widely varying runtime speeds, especially the first time it is run.

In [1]:

from datetime import datetime, timedelta
import os
import glob
import shutil

import param
import paramnb
import numpy as np
import xarray as xr
import geoviews as gv
import holoviews as hv
import quest
import earthsim.gssha as esgssha
import earthsim.gssha.model as models
import cartopy.crs as ccrs

from earthsim.gssha import download_data, get_file_from_quest
from holoviews.streams import PolyEdit, BoxEdit, PointDraw, CDSStream
from holoviews.operation.datashader import regrid, shade
from earthsim.io import save_shapefile, open_gssha, get_ccrs

regrid.aggregator = 'max'

hv.extension('bokeh')
%output holomap='scrubber' fps=2

In [2]:

shutil.rmtree('./vicksburg_south',ignore_errors=True)

Configure model parameters ¶

In [3]:

model_creator = esgssha.CreateGSSHAModel(name='Vicksburg South Model Creator',
                                        mask_shapefile='../data/vicksburg_watershed/watershed_boundary.shp',
                                        grid_cell_size=90)
paramnb.Widgets(model_creator,initializer=paramnb.JSONInit())

Draw bounds to compute watershed ¶

Allows drawing a bounding box and adding points to serve as input to compute a watershed:

In [4]:

%%opts Polygons [width=900 height=500] (fill_alpha=0 line_color='black')
%%opts Points (size=10 color='red')
tiles = gv.WMTS('http://c.tile.openstreetmap.org/{Z}/{X}/{Y}.png',
                crs=ccrs.PlateCarree(), extents=(-91, 32.2, -90.8, 32.4))
box_poly = hv.Polygons([])
points = hv.Points([])
box_stream = BoxEdit(source=box_poly)
point_stream = PointDraw(source=points)
tiles * box_poly * points

Out[4]:

In [5]:

if box_stream.element:
    element = gv.operation.project(box_stream.element, projection=ccrs.PlateCarree())
    xs, ys = element.array().T
    bounds = (xs[0], ys[1], xs[2], ys[0])
    print("BOUNDS", bounds)
    
if point_stream.element:
    projected = gv.operation.project(point_stream.element, projection=ccrs.PlateCarree())
    print("COORDINATE:", projected.iloc[0]['x'][0], projected.iloc[0]['y'][0])

Inspect and edit shapefile ¶

The plot below allows editing the shapefile using a set of tools. The controls for editing are as follows:

Double-clicking the polygon displays the vertices
After double-clicking the point tool is selected and vertices can be dragged around
By tapping on a vertex it can be selected, tapping in a new location while a single point is selected inserts a new vertex
Multiple points can be selected by holding shift and then tapping or using the box_select tool
Once multiple vertices are selected they can be deleted by selecting the point editing tool and pressing backspace

In [6]:

%%opts Shape [width=900 height=500 tools=['box_select']] (alpha=0.5)
mask_shape = gv.Shape.from_shapefile(model_creator.mask_shapefile).last
tiles = gv.WMTS('http://c.tile.openstreetmap.org/{Z}/{X}/{Y}.png')
vertex_stream = PolyEdit(source=mask_shape)
tiles * mask_shape

Out[6]:

If any edits were made to the polygon in the plot above we save the watershed_boundary.shp back out and redisplay it to confirm our edits were applied correctly:

In [7]:

%%opts Shape [width=600 height=400] (alpha=0.5)
if vertex_stream.data:
    edited_shape_fname = './vicksburg_watershed_edited/watershed_boundary.shp'
    dir_name = os.path.dirname(edited_shape_fname)
    if not os.path.isdir(dir_name): os.makedirs(dir_name)
    save_shapefile(vertex_stream.data, edited_shape_fname, model_creator.mask_shapefile)
    model_creator.mask_shapefile = edited_shape_fname
    mask_shape = gv.Shape.from_shapefile(edited_shape_fname).last
mask_shape = mask_shape.opts() # Clear options
mask_shape

Out[7]:

Configure simulation parameters ¶

In [8]:

sim = esgssha.Simulation(name='Vicksburg South Simulation', simulation_duration=60*60,
                          rain_duration=30*60, model_creator=model_creator)
paramnb.Widgets(sim,initializer=paramnb.JSONInit())

Create the model ¶

Note that the above code demonstrates how to collect user input, but it has not yet been connected to the remaining workflow, which uses code-based specification for the parameters.

In [9]:

if sim.model_creator.project_name not in quest.api.get_collections():
    quest.api.new_collection(sim.model_creator.project_name)

In [10]:

paramnb.Widgets(sim.model_creator,initializer=paramnb.JSONInit())

In [11]:

# temporary workaround until workflow cleanup/parameterization is done
if sim.model_creator.project_name == 'test_philippines_small':
    sim.model_creator.roughness = models.GriddedRoughnessTable(
        land_use_grid=get_file_from_quest(sim.model_creator.project_name, sim.land_use_service, 'landuse', sim.model_creator.mask_shapefile),
        land_use_to_roughness_table='./philippines_small/land_cover_glcf_modis.txt')
else:    
    sim.model_creator.roughness = models.GriddedRoughnessID(
        land_use_grid=get_file_from_quest(sim.model_creator.project_name, sim.land_use_service, 'landuse', sim.model_creator.mask_shapefile),
        land_use_grid_id=sim.land_use_grid_id)

sim.model_creator.elevation_grid_path = get_file_from_quest(sim.model_creator.project_name, sim.elevation_service, 'elevation', sim.model_creator.mask_shapefile)

INFO:rasterio._base:Failed to auto identify EPSG: 7

... ... .tif format raster saved at /home/travis/.local/share/quest/projects/default/vicksburg_south/NLCD2011_LC_N30W090.tif
retrieving raster availability from https://www.sciencebase.gov/catalog/items?parentId=4f70aa9fe4b058caae3f8de5&filter=tags%3DIMG&max=1000&fields=webLinks%2Cspatial%2Ctitle&format=json&filter=spatialQuery%3D%7Bwkt%3A%22POLYGON+%28%28-180.0+90.0%2C-180.0+-90.0%2C180.0+-90.0%2C180.0+90.0%2C-180.0+90.0%29%29%22%2Crelation%3A%22intersects%22%7D
retrieving raster availability from https://www.sciencebase.gov/catalog/items?filter=tags%3DIMG&filter=spatialQuery%3D%7Bwkt%3A%22POLYGON+%28%28-180.0+90.0%2C-180.0+-90.0%2C180.0+-90.0%2C180.0+90.0%2C-180.0+90.0%29%29%22%2Crelation%3A%22intersects%22%7D&max=1000&format=json&fields=webLinks%2Cspatial%2Ctitle&parentId=4f70aa9fe4b058caae3f8de5&offset=1000
... ... .img format raster saved at /home/travis/.local/share/quest/projects/default/vicksburg_south/USGS_NED_13_n33w092_IMG.img
... ... .img format raster saved at /home/travis/.local/share/quest/projects/default/vicksburg_south/USGS_NED_13_n33w091_IMG.img

In [12]:

model = sim.model_creator()

In [13]:

# add card for max depth
model.project_manager.setCard('FLOOD_GRID',
                              '{0}.fgd'.format(sim.model_creator.project_name),
                              add_quotes=True)
# Add time-based depth grids to simulation
"""
See: http://www.gsshawiki.com/Project_File:Output_Files_%E2%80%93_Required

Filename or folder to output MAP_TYPE maps of overland flow depth (m) 
every MAP_FREQ minutes. If MAP_TYPE=0, then [value] is a folder name 
and output files are called "value\depth.####.asc" **
"""

model.project_manager.setCard('DEPTH', '.', add_quotes=True)
model.project_manager.setCard('MAP_FREQ', '1')

# add event for simulation (optional)
"""
model.set_event(simulation_start=sim.simulation_start,
                simulation_duration=timedelta(seconds=sim.simulation_duration),
                rain_intensity=sim.rain_intensity,
                rain_duration=timedelta(seconds=sim.rain_duration))
"""
# write to disk
model.write()

Review model inputs ¶

Load inputs to the simulation ¶

In [14]:

name = sim.model_creator.project_name
CRS = get_ccrs(os.path.join(name, name+'_prj.pro'))

roughness_arr = open_gssha(os.path.join(name,'roughness.idx'))
msk_arr = open_gssha(os.path.join(name, name+'.msk'))
ele_arr = open_gssha(os.path.join(name, name+'.ele'))

roughness = gv.Image(roughness_arr, crs=CRS, label='roughness.idx')
mask = gv.Image(msk_arr, crs=CRS, label='vicksburg_south.msk')
ele  = gv.Image(ele_arr, crs=CRS, label='vicksburg_south.ele')

Shapefile vs. Mask ¶

In [15]:

tiles * regrid(mask) * mask_shape

Out[15]:

Elevation ¶

In [16]:

tiles * regrid(ele) * mask_shape

Out[16]:

Roughness ¶

In [17]:

tiles * regrid(roughness) * mask_shape

Out[17]:

Run Simulation ¶

In [18]:

from gsshapy.modeling import GSSHAFramework

In [19]:

# TODO: how does the info here relate to that set earlier?

# TODO: understand comment below
# assuming notebook is run from examples folder
project_path = os.path.join(sim.model_creator.project_base_directory, sim.model_creator.project_name)
gr = GSSHAFramework("gssha",
                    project_path,
                    "{0}.prj".format(sim.model_creator.project_name),
                    gssha_simulation_start=sim.simulation_start,
                    gssha_simulation_duration=timedelta(seconds=sim.simulation_duration),
                    # load_simulation_datetime=True,  # use this if already set datetime params in project file
                   )

# http://www.gsshawiki.com/Model_Construction:Defining_a_uniform_precipitation_event
gr.event_manager.add_uniform_precip_event(sim.rain_intensity, 
                                          timedelta(seconds=sim.rain_duration))

gssha_event_directory = gr.run()

Visualizing the outputs ¶

Load and visualize depths over time ¶

In [20]:

depth_nc = os.path.join(gssha_event_directory, 'depths.nc')
if not os.path.isfile(depth_nc):
    # Load depth data files
    depth_map = hv.HoloMap(kdims=['Minute'])
    for fname in glob.glob(os.path.join(gssha_event_directory, 'depth.*.asc')):
        depth_arr = open_gssha(fname)
        minute = int(fname.split('.')[-2])
        # NOTE: Due to precision issues not all empty cells match the NaN value properly, fix later
        depth_arr.data[depth_arr.data==depth_arr.data[0,0]] = np.NaN
        depth_map[minute] = hv.Image(depth_arr)

    # Convert data to an xarray and save as NetCDF
    arrays = []
    for minute, img in depth_map.items():
        ds = hv.Dataset(img)
        arr = ds.data.z.assign_coords(minute=minute)
        arrays.append(arr)
    depths = xr.concat(arrays, 'minute')
    depths.to_netcdf(depth_nc)
else:
    depths = xr.open_dataset(depth_nc)

depth_ds = hv.Dataset(depths)
depth_ds.data

Out[20]:

<xarray.Dataset>
Dimensions:  (minute: 60, x: 464, y: 400)
Coordinates:
  * x        (x) float64 6.67e+05 6.671e+05 6.672e+05 6.673e+05 6.674e+05 ...
  * y        (y) float64 3.548e+06 3.548e+06 3.548e+06 3.548e+06 3.548e+06 ...
  * minute   (minute) int64 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 ...
Data variables:
    z        (minute, y, x) float64 nan nan nan nan nan nan nan nan nan nan ...

Now that we have a Dataset of depths we can convert it to a series of Images.

In [21]:

%%opts Image [width=600 height=400 logz=True xaxis=None yaxis=None] (cmap='viridis') Histogram {+framewise}
regrid(depth_ds.to(hv.Image, ['x', 'y'])).redim.range(z=(0, 0.04)).hist(bin_range=(0, 0.04))

Out[21]:

We can also lay out the plots over time to allow for easier comparison.

In [22]:

%%opts Image [width=300 height=300 logz=True xaxis=None yaxis=None] (cmap='viridis')
regrid(depth_ds.select(minute=range(10, 70, 10)).to(hv.Image, ['x', 'y']).redim.range(z=(0, 0.04))).layout().cols(3)

Out[22]:

Flood Grid Depth ¶

(Maximum flood depth over the course of the simulation)

In [23]:

%%opts Image [width=600 height=400] (cmap='viridis')
fgd_arr = open_gssha(os.path.join(gssha_event_directory,'{0}.fgd'.format(sim.model_creator.project_name)))
fgd = gv.Image(fgd_arr, crs=CRS, label='vicksburg_south.fgd').redim.range(z=(0, 0.04))
regrid(fgd, streams=[hv.streams.RangeXY]).redim.range(z=(0, 0.04))

Out[23]:

Analyzing the simulation speed ¶

In [24]:

%%opts Spikes [width=600]
times = np.array([os.path.getmtime(f) for f in glob.glob(os.path.join(gssha_event_directory, 'depth*.asc'))] )
minutes = (times-times[0])/60
hv.Spikes(minutes, kdims=['Real Time (minutes)'], label='Time elapsed for each minute of simulation time') +\
hv.Curve(np.diff(minutes), kdims=['Simulation Time (min)'], vdims=[('runtime', 'Runtime per minute simulation time')]).redim.range(runtime=(0, None))

Out[24]:

Here if the "spikes" are regularly spaced, simulation time is regularly scaled with real time, and so you should be able read out the approximate time to expect per unit of simulation time.

GSSHA Workflow ¶